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Abstract. - We investigate the modifications brought about by the linear connectivity among 
charges in the classical Thomson problem. Instead of packing with local hexagonal order inter- 
sperced with topological defects, we find charge distributions with helical symmetry wound around 
the surface of the sphere. This finding should have repercussions in the viral packing and macroion 
adsorption theories. 



Introduction. - Electrostatics can be very unintuitive and can lead to quite unexpected 
features in a variety of physical, chemical and biological contexts [1]. The problem of packing 
charges on a sphere, generally referred to as the Thomson problem [2], has recently received a 
lot of attention and has been analysed via detailed direct numerical evaluation [3] as well as via 
approximate analytical theories [4] . What is clear is that for large enough number of charges the 
smallest electrostatic energy configurations have packing with local hexagonal order intersperced 
with topological defects [4]. A possible question here is: what happens if one linearly connects 
the charges and thus considers a Hnearly connected Thomson problem. 

The linearly connected Thomson problem should be directly relevant for packing or adsorption 
of polyelectrolytes, i.e. charged polymers, on or to spherical surfaces. Adsorption of polyelec- 
trolytes to spheres has been extensively simulated in the context of polyelectrolyte - macroion 
interactions in colloid physics [5]. Various configurations of adsorbed polyelectrolytes are seen, 
depending on the interaction parameters and the stiffness of the chain. Though they all cor- 
respond in one sense or another to a free energy minimisation as in the case of histone - DNA 
interactions [6], the phase space is simply too big in order to attempt some kind of systematic 
classification of shapes of the adsorbed polyelectrolyte chains. An "inverted" problem to the one 
just described is represented by the investigations of the nature of packing of DNA inside small 
bacteriophage capsids that has recently received a lot of attention since the nature of packing of 
DNA inside the capsid determines the ejection force that can be measured directly [7]. The only 
general and consistent conclusion stemming from these studies seems to be that DNA close to 
the capsid wall is more ordered than deep within the capsid. Recent cryomicroscopy on epsilonl5 
bacteriophage [8] gives ample experimental support for this type of order gradient within the 
viral capsid. At small packing fractions, where most or indeed all of DNA is expelled towards 
the inner surface of the capsid, the packing problem is directly related to the linearly connected 
Thomson problem. 
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Fig. 1 - This Figure illustrates the geometry used in this paper. Three charges labelled A, B and C lie 
on a surface of the sphere, whose centre is at Q. Point P is such that AP and PC are perpendicular to 
the line QB. The angles Oi used in this work correspond to tt — {APC. 



In order to shed some light onto packing of connected objects into confined spaces with 
spherical topology, we performed zero-temperature numerical investigations of the connected- 
charges Thomson effect. Apart from the fundamental significance of this problem [2] it could 
also be seen as a model of polyelectrolyte packing either on the outside or the inside of a sphere 
in the cases where the entropy of the chain can be neglected. 

Model. - Our model system is composed of N+1 elementary charges connected by N rigid 
bonds of length ^ on a unit 3D sphere. Each rod subtends and angle 4> = 2sin^^(i'/2) with the 
centre of the sphere. The position of the charges are given by + 1 vectors Xj (z = 1...N + 1), 
satisfying |xi| = 1 and |xj— Xj_i| = £. These constraints reduce the number of degrees of freedom 
to iV + 2. Imposing rotational invariance further restricts the system to just — 1 degrees of 
freedom. The total length of the chain is L = N(l). The Coulomb energy of this system is given 

by 

i,j>i 

It is useful for comparison to divide this value by the minimum energy the chain would have if 
it was not constrained to lie on a sphere. In that case, the minimum energy configuration would 
be a completely straight chain, whose energy is given by 

^ 1 1 

i,j>i ' 

where ^ is the standard digamma function and 7 is the Euler's constant. 

In the limit of large A^, the minimal energy of an unconstrained Coulomb chain reduces to 

Eo{L)r^jN\ogN (3) 

This has a fairly obvious explanation: for a long enough chain, most of the potential energy 
comes from the nearest neighbour interactions and therefore scales as N, while the Coulomb 
interactions between more remote charges contribute only logarithmically in N . This result is a 
signature of a completely extended Coulomb chain. 
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Fig. 2 - This figure shows the reduced potential energy as a function of the total chain length L for 
the four cases considered here: 4> — 7v/4 (top), vr/lO, 7r/20, -k/40 (bottom). Solid line corresponds to the 
annealing in terms of ai, while the dotted thick line corresponds to annealing in terms of l3i. For large 
values of A'', the computations become prohibitively expensive. 



We define the reduced potential energy by 

PiL) = 



EjL) 



(4) 



which measures the ratio between the Coulomb self-energy of a constrained and completely 
extended chain. 

A natural description of this system is by — 1 angles ji subtended between rigid rods linking 
i — 1th, zth and i + 1th charge, i.e. cos7i — (x^ — Xi_i) • (x^+i — x^). However, by projecting the 
mid-point down to the plane perpendicular to x^ and going through Xi_i and x^+i one arrives 
to a mathematically more convenient description in terms of the angles at . We will define N —I 
angles Ui as 



COStti = (c0S(/)Xi — Xi_i) • (Xi+i — COS0X 

The geometry of this choice of variables is illustrated on Figure 
Using straightforward geometry one can show that 

cos + fi COS ai + Pi sin ai , 



with 



and 



fi = Xj cos (j) - Xj 



X,,;Xf,; 



(5) 

(6) 
(7) 



(8) 

By construction x^+i lies on a sphere if both Xi_i and x^ lie on a sphere. An = configuration 
corresponds to a linear, completely extended chain, while ai = const. ^ configuration results 
in charges being distributed along an arc on a spherical surface. 

We will also use — 1 coefficients Pi that can be derived simply from a polynomial expansion 
of ai around the midpoint of the chain and are defined by 

j = l...JV-l ^ ' ^ 
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Fig. 3 - Typical minimum energy configurations for charges on a unit sphere. The l.h.s. corresponds to 
a typical solutions obtained with ai minimisation, while the r.h.s. is one obtained by Pi minimisation. 
Both are very similar (within better than 0.1%) in terms of the actual energy. 

Description of the chain by Pi or Ui is only a matter of convenience. 

Saff and Kuijlaars [10] have shown that the "general spiral set" can be used to generate 
uniformly distributed points on a sphere. Their formalism can not however be directly applied 
to our problem, since the fixed Hnks constrain the system further. Though we are not concerned 
explicitly with uniform distributions of points on a sphere, one can, nevertheless, reasonably 
expect that the minimal solution for the Thomson problem in the constrained case will exhibit 
helical behaviour. This is a hypothesis that we will try to investigate. 

In order to set the terms of discussion, we recall the equation of a helix on a cylinder as 
given hy z = ct and cf) — t. This can be generalised to a helix on a surface of a sphere by either 
identifying z with cos 9 or 9. We will denote the first case as embedded helix: 

e = cos"^(ci) 

= t (10) 

and the second one as spherical helix: 

9 = ct 

4> = t, (11) 

where 9 and <f) are the ordinary spherical coordinates. 

Results. - We have programmed software that finds the minimum energy configuration 
using two methods, working in either ai space or Pi space. 

In the first case we anneal by angles ai using a version of genetic algorithm used in [9] for a 
classical Thomson problem. We start with a population of P "parents". We create P{P — l)/2 
"off-springs" by taking first N/2 Ui values from the first parent and the remaining a^s from the 
second parent for all possible combinations of parents. We then use a simple numerical minimiser 
(downhill simplex) to find the local potential energy minimum corresponding to each parent. The 
best P solutions are kept and the process is repeated until a satisfactory convergence is reached. 
We typically work with 12 initial parents and 8 generations. 
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Fig. 4 - This figure shows the difference in reduced potential energy P between solutions found by full 
numerical minimisation and those found by constraining charges to lie on either embedded helix (solid) 
or spherical helix (dashed). Lines were offset vertically by i x 0.1 and correspond the four cases considered 
here: = n/4 (top), tt/IO, tt/20, tt/AO (bottom). The dotted lines represent the zero difference value. 
For large values of A'', the computations become prohibitively expensive. 



In the second case we find the best solution by minimising the energy expressed in terms of 
parameters (3i. Genetic algorithm is not useful here as each Pi affects the entire chain and so 
good local features cannot be preserved. However, it is an extremely good space for optimisation 
and thus finds a good solution very quickly. After trying several options we have settled with an 
annealing algorithm that is a series of linear minimisers in randomly chosen parameters where 
the solution is assumed to be bracketed around the last optimal point and the width of the 
bracket is decaying exponentially until a desired convergence is reached. We typically try ~ 200 
different starting positions and then choose the best. 

Finally, we have also minimised the energy of the system, whose charges are additionally 
constrained to He on either embedded or spherical helix. Since this is a one-parameter model, 
the straight numerical minimiser performs satisfactorily. The purpose of this exercise was to 
check how well the true minimal configuration can be approximated by either helix. 

We have run the minimiser with four values of (fi: 7r/4, tt/IO, 7r/20 and 7r/40 and for values 
of L ranging between 7r/2 and IOtt. 

In Figure|2lwe plot P{L) as a function of the total chain length. It is immediately clear that 
both annealing procedures lead to stable solutions that are very similar in terms of energies for 
large enough values of the length of the chain. The well discernible minimum in the P{L) is easy 
to understand as follows: the energy Eq{L) is dominated by the nearest-neighbour contribution 
to the overall potential. The penalty payed by constraining the chain to He on a sphere comes 
primarily from the "bending of the chain". However, when L > it, the E{L) receives another 
contribution from interaction between charges that lie far apart on the chain but a forced to lie 
close in physical coordinates due to wrapping of the chain around the sphere. Therefore, one can 
expect a minimum in P at around L/tt ~ 1 and this is exactly what we observe. However, we 
note that for a long enough chain, the reduced potential energy begins to scale nearly linearly 
with the chain length. 

Naively, one would expect that the energy of the electrons on the chain wrapped on a sphere 
would scale with N roughly as N'^ and so the reduced potential energy P should scale as N/ log N 
for large N. Our data do not extend high enough in TV to discriminate unequivocally between 
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linear and Nj \ogN scaling, though the latter indeed provides a marginally better fit to the data. 
Interestingly enough, this is the true even for cases when values of Eg have not yet reached the 
iVlogiV scaling. 

We also note, that although very similar in energy, the solutions found by the Pi anneahng 
procedure are always helical, while the solutions in terms of the ai annealing procedure are, for 
large enough N, configurations that are locally (on a sphere, not along the chain) helical but are 
irregular on large scales. This is illustrated in the Fig. What we are likely seeing is a plethora 
of metastable configurations with nearly minimal energy, similar to those found for a classical 
Thomson problem. This is supported by two observations: firstly, making minimiser convergence 
criteria stronger does not improve these solutions, implying that they are indeed meta-stable 
and secondly, all final generation solutions in the genetic algorithm have very similar energies 
(to within a few pro-miles) and are different on runs with different random seeds, implying 
that the family of those solutions is quite large. These nearly minimal configurations are never 
found by the (3i annealing as this would require a coherent tuning of several Pi coefficients to 
achieve a preferred local behaviour. The existence of a large number of these nearly minimal 
configurations also implies that unless a specific mechanism is in place, the chain is very likely 
to land in a complex nearly-perfect meta-stable state rather than the helical structure with 
long-scale ordering. 

Next let us discuss the solutions which come from constraining the charges to He on either 
of the two helices defined in Eqs. ifTTll and itTTHl . Results are plotted in the Figure 21 It can be 
seen that for chains with 4> ^ tt, the spherical helix is an increasingly better approximation to 
the true minimum for large enough L. It thus appears that we are able to construct explicitly 
those configurations that are close to optimal. This confirms our expectations based on the 
"generalised spiral set" of Saff and Kuijlaars [10] that solutions with heHcal symmetry should 
describe minimal energy states for the connected-charges Thomson problem. 

Conclusions. - In this paper we have analysed the connected-charges Thomson problem by 
brute-force numerical minimisation of the potential energy. The problem is completely specified 
by two parameters: L, the total length of the chain, and = L/N, that defines the inter-charge 
separation. Note that the links between the charges are inextensible. 

We have found that for all cases considered, the reduced potential energy has a minimum at 
L/tt ^ 1 and scales as N or N log N (the range of iVs investigated does not extend to large enough 
N so that one could differentiate between the two scaling forms) . This implies furthermore that 
the potential energy of the chain scales as for large N as expected from standard arguments. 
However, the relation is true even for N that are so low that Ef){N) has not reached the A^logiV 
scaHng yet. 

We have heuristically attempted to model the minimum energy configuration with a helix 
on a sphere. We have found that the cylindrical helix projected onto a sphere by a simple 
transformation z — > 6*, thus creating a spherical helix, leads to charge configurations that seem 
to be very close to the numerically obtained minimum for L > tt and 0^1. The solution is 
likely to continue in such manner until the distance between points that are far apart along the 
helix but close in real space becomes important. This would result in a frustrated helix that 
might hide new and exciting phenomenology. Unfortunately, the presently available computers 
do not reach the computational power necessary to answer this question, at least not with the 
methods provided in this paper. 

What is the physical relevance of our results? Adsorption of fiexible polyelectrolytes onto 
spherical macroions (colloidal interactions) and configurations of polyelectrolytes within spherical 
cavities (viral packing) both belong to the type of problems which have been analysed in this 
paper on their bare bones level. Though in general the configurations of real polyelectrolyte chains 
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depone! on many parameters and are difRciilt to classify [11], our analysis provides an asymptotic 
limiting configuration that should be discernible in the case where adsorption interaction onto a 
spherical macroion is strong enough, or the repulsive interactions along the chain in a spherical 
cavity are sufficiently long-ranged. An especially beautiful illustration of our theory is provided 
by recent elucidation of the DNA packing inside the epsilonl5 capsid [8] that shows tight helical 
winding on the outer boundary and the progressive loosening of the orientational order as one 
moves away from the boundary towards the interior of the capsid. 
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